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Abstract 

Simulations of biophysical systems inevitably include steps that correspond to time 
integrations of ordinary differential equations. These equations are often related 
to enzyme action in the synthesis and destruction of molecular species, and in the 
regulation of transport of molecules into and out of the cell or cellular compartments. 
Enzyme action is almost invariably modeled with the quasi-steady-state Michaelis- 
Menten formula or its close relative, the Hill formula: this description leads to 
systems of equations that may be stiff and hard to integrate, and poses unusual 
computational challenges in simulations where a smooth evolution is interrupted by 
the discrete events that mark the cells' lives. This is the case of a numerical model 
(Virtual Biophysics Lab - VBL) that we are developing to simulate the growth of 
three-dimensional tumor cell aggregates (spheroids). The program must be robust 
and stable, and must be able to accept frequent changes in the underlying theoretical 
model: here we study the applicability of known integration methods to this unusual 
context and we describe the results of numerical tests in situations similar to those 
found in actual simulations. 
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1 Introduction 



Most mathematical models of biochemical and biophysical processes in cells 
are described by nonlinear differential systems that cannot be handled with 
analytical methods and require numerical integrations [lll2] . However the high 
speed of computers and the sophisticated computational methods that are 
available today are powerful tools that allow the numerical exploration of these 
exceedingly interesting dynamical systems. This also suggests that eventually 
the biophysical models will no longer be analytic, but mostly computational; 
and indeed we are now developing a program that simulates tumor spheroids 
(VBL, Virtual Biophysics Lab) [3llif5|6] . which includes a reduced but still 
quite complex description of the biochemistry of individual cells, plus many 
diffusion processes that bring oxygen and nutrients into cells and metabolites 
into the environment. 

In a numerical model like VBL, each cell is described by a reduced metabolic 
network, and by other mechanisms that include both discrete deterministic 
and stochastic events. The description is thus mixed, with smooth evolutions 
interspersed with discrete steps. The exchange of molecules with the surround- 
ing environment means that transport into and out of cells is closely linked 
with diffusion processes that involve the whole cluster of cells, and finally lead 
to a very large set of (time) differential equations. Simulation steps between 
discrete events require the integration of nonlinear differential equations that 
describe the individual cell's clockwork, and the integration of the diffusion 
equations. These integrations are carried out under widely different condi- 
tions, in a changing environment, and for this reason they need algorithms 
that are both unconditionally stable and free from unwanted artifacts. These 
conditions are not always fulfilled in the existing literature, and we feel that 
a detailed understanding of the underlying mathematical and computational 
bases may be important not just for us, but for other workers in the field of 
systems biology as well. 

Eventually the whole simulation program must run smoothly, and it must be 
free of stability problems. Thus it is very important to ensure that step by 
step integrations are stable and do not bring the biochemical variables into 
unphysical regions (e.g., no concentration must ever become negative). 



2 Biophysical models 

Biophysical models are not arbitrary dynamical systems: indeed the phenom- 
ena of life are characterized by remarkable stability properties, as most of 
them display either homeostasis or (nearly) stable limit cycles. At the basic 
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reaction level, many processes obey a straightforward enzyme kinetics, regu- 
lated by the well-known Michaelis-Menten equation pQ, or by a variant that 
applies to cooperative processes, the Hill equation [7]. 

Here we start with a special case that displays generic features shared by 
many other processes, namely the transport of glucose into and out of cells 
and its conversion into glucose-6-phosphate (G6P; this is part of the reduced 
metabolic model incorporated in our simulator of cell metabolism, growth and 
proliferation [if5] ). We begin with the equations for a single cell in a stable 
environment: 
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Here the dynamical variables are the glucose concentration inside the cell [Gm] , 
and the concentration in the extracellular space [Gextra]', the environmental 
glucose concentration [Genv] is fixed and is one of the model parameters. The 
introduction of the extracellular space may look like a useless nuisance in this 
case, but it is justified on two different, and both important, grounds. The ex- 
tracellular space is absolutely needed in multicellular systems, because in that 
case it becomes the scaffolding that allows the diffusion of many substances 
in closely bound cell clusters (and corresponds to an actual biological entity, 
the free space available in the extracellular matrix piQfTO] ). and a consistent 
treatment requires its extension also towards the environment; secondly, the 
extracellular space, even in this simple case, acts as the cell-environment in- 
terface, and corresponds to the buffer region actually observed in diffusion 
measurements in tumor spheroids [llj. 

We assume a spherical cell of radius r ^ 5yum, so that the surface area is 
S = Airr'^ pa 3 ■ lO'^^m^, and the volume is Vin = (4/3)7rr^ ^ 5 • lO'^'^m^. 
The first two terms on the rhs of equation ([T]) correspond to the transport of 
glucose from the extracellular space into the cell and the reverse process of 
transport from the inside of the cell back to the extracellular space (enacted 
by the reversible GLUT transporters, see [T2IT3] and references therein); since 
this transport process is facilitated P^PH] we describe it with two Michaelis- 
Menten terms, where the maximum transport speed fmax,i is proportional 
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to the cell's surface ji|12m3j . and Ki is an experimental parameter [Ij.The 
other two terms correspond to a cascade of two Michaelis-Menten processes 
that involve the enzymes glucokinase and hexokinase, and is regulated by 
parameters that do not depend on cell size [H]. 

The second equation (eq. (|2])), in addition to the GLUT- mediated transport 
terms, includes a standard diffusion term that corresponds to the diffusion of 
glucose from the surrounding environment into the extracellular space. This 
term is a discretization of glucose flux across the surface S (which is approxi- 
mately the area of the cell-extracellular space interface); the distance used to 
compute the concentration gradient is the thickness of the extracellular layer, 
A ^ 0.2/im. 

Although equations like ([T]) and (|2| are usually cast in terms of concentrations, 
we prefer to use masses, as this makes total mass conservation stand out 
clearly: 
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We postpone the listing of model parameters to section |4| where we summarize 
the results of numerical tests. 



3 Integration methods 



There is a broad and comprehensive literature on integration methods, and 
at first it may seem that in large-scale biophysical simulations the choice is 
only limited by the required accuracy and by the available computing power, 
however it is not so. In a simulation like VBL [IfS] there are many cells and 
the number of equations is quite large: our aim is to simulate at least one 
million cells, which corresponds to a spherical cell cluster with a diameter of 
one millimeter. In VBL we make several simplifying assumptions and one of 
them is that each cell has its own extracellular space: this means that in such 
a system - for each molecular species which diffuses in the cell cluster - there 
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are about two million local concentration variables (for each cell there is an 
intracellular and an extracellular concentration) and a corresponding number 
of equations. 

The simulation includes random events as well (such as the partitioning of 
organelles at the time of cell division), so that numerical inaccuracies are ab- 
sorbed in this much larger randomness and turn out to be much less important 
than they are in strictly deterministic systems. 

Moreover, the life of cells is marked by discrete events: this discreteness pro- 
duces endless transient responses in the deterministic part of the simulation 
algorithm, and the most evident discrete process - cell division - also changes 
the number of variables and equations. 

Last but not least, there is a wide range of volumes (cells and extracellular 
spaces), and this leads to very stiff systems of equations: in fact while the 
metabolic processes can be quite slow (with characteristic times of the order 
of thousands of seconds), a single extracellular space has a very small volume 
and the time it takes to fill or empty this volume is very short (of the order of 
10-100 /is). Thus the characteristic times span eight orders of magnitude, and 
the system is quite stiff. This means that any adopted integration method must 
be very stable, and that it should be able to overlook very short transients 
that do not mean much from a biophysical point of view, and it should also 
accept long time steps without crashing or producing unphysical results. 

All this proves quite challenging, as it poses conflicting requests on the in- 
tegration algorithms. The usual considerations on stability [T5]ll6l)17j suggest 
that implicit methods should perform better than explicit methods, and this 
also means that a robust algorithm must sacrifice speed, as implicit methods 
require several function evaluations. 

Multistep methods p^l7|ll8] may appear to be a smart choice, because they 
can incorporate information from previous simulation steps and thus reduce 
the number of required function evaluations: however numerical tests per- 
formed in anticipation of the present work have shown that explicit Adams- 
Bashforth multistep methods are highly unstable in this context, and that the 
fourth-order Adams-Moulton implicit multistep method does not fare much 
better as it seems to be very sensitive to the quality of the algorithm initial- 
ization steps (that must necessarily be performed with another integration 
algorithm), and unfortunately the discrete events lead to a constant flow of 
equation initialization processes. 

Finally we have decided to include the following methods in the tests summa- 
rized in the next section: 

• the default integrator of Mathematica [12] , as a benchmark algorithm; 
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• the implicit Euler method; 

• the imphcit trapezoidal method (which corresponds to the second-order 
Adams-Moulton algorithm) . 

Although it cannot be included as it is in a simulation program, we have de- 
cided to choose the default setting of the Mathematica instruction NDSolve 
as a benchmark since it incorporates a well-tuned mixture of standard algo- 
rithms not included in our very short list [19], and thus provides a qualitative 
guide to possible better choices. Here simphcity is an important bonus, since 
the integration algorithm must be incorporated in a larger structure and the 
biophysical model may change frequently, as new biochemical paths are incor- 
porated: however we have not included the simple explicit Euler method in 
the list, because it is well-known to fail badly in the case of stiff equations, 
unless the step size is exceedingly small (and yet this method is widely used 
in many similar contexts, see, e.g, [2D])- The many important Runge-Kutta 
(RK) methods are also missing, because the stable implicit RK methods are 
not easily adopted for inclusion in a simulation program like VBL (or other 
similar numerical stepping schemes), and also because the sophisticated inte- 
grator of Mathematica [19j actually includes a choice of RK methods - both 
explicit and implicit - and its performance may suggest corrections to our 
selection of algorithms. 

Since stability - and not precision - is the main feature that an integrator 
must have in this context, it would be important to set the general problem 
in a form such that the standard theorems on stability can be applied [16] . 
Unfortunately it is impossible to specify exactly the operating conditions of 
the simulator, especially because the number of equations is variable - as 
cells proliferate - and because the theoretical model must be open to changes 
as our understanding of the biophysical and biochemical processes evolves. 
However there are a couple of general considerations that can help: the first 
is that the systems of equations must be autonomous, time cannot appear 
explicitly. The second is that the equations describe the evolution of masses 
and concentrations, i.e. quantities that must be non-negative, and this means 
that the generic structure of the equations for non-negative quantities must 
be 

— = A{xi, ...,XN]t)- B{xi, ...,xn; t)xk (5) 



where A and B are non-negative functions that express respectively production 
and consumption/ destruction of the quantity Xk- The consumption/destruction 
term must be proportional to Xk to ensure the continued non-negativity of the 
solution (if there were not such a proportionality, then the derivative could 
be nonvanishing and negative for values of Xk close to zero, and thus lead to 
negative Xk)- Indeed the proportionality is often related to Michaelis-Menten 



6 



terms with a generic structure Xk/ {Km + Xk)- A differential system such as ([s]) 
can be locally approximated by the simpler equations 



dxk 

~dt 



Aq - B^Xk 



(6) 



which shows that A-stable integration methods should normally be sufficient 
in this complex context [16]. 



4 Numerical results on the test model of section [2] 

We have chosen parameter values derived from experiment or from fits of 
experimental data and they are listed in table 1; all parameters are further 
explained and referenced in [Hl5] . except D |21] and A. The approximate value 
of A is derived from the estimate that the extracellular matrix amounts to 
about 20% of the total volume of human tissue [10] and assuming a free space 
fraction of about 50% [8j (we note that this may actually be an overestimate): 
then we obtain A 0.5 {Q.2hVin) /S ^ 0.2/im. 

Using the values in the table we can find the equilibrium values that correspond 
to drriin/dt = and drriextra/dt = 0, i.e., [Gin] ~ 0.043 kg m"^ and [Gextra] ~ 
0.90 kg m~^: this means that because of conversion into G6P, the glucose level 
in this isolated cell is considerably lower than in the surrounding environment. 

As explained above, the system of equations ^ and ^ is very stiff, and this 
can be easily visualized setting an off-equilibrium initial glucose concentration 
inside the cell. Some of the glucose is quickly converted into G6P, but part 
of it either seeps out or enters the cell, and changes the concentration in the 
extracellular space. Since the extracellular volume is very small it reacts very 
quickly to this sudden change, and exchanges glucose with the environment, 
exhibiting a sharp transient. The approximate duration of this transient can 
be estimated expanding equation Q for very low values of the masses: since 



with the parameters of table 1, this gives niextrait) ~ ^ + -Bexp (— t/r), with 
T 60/is. It is also important to note that this very short time is entirely due to 
the diffusion term in equation ([t]): without the diffusion term the characteristic 
time r would grow to more than 27 s. 
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parameter 


value 




r 


• iO m 


cell radius 


D = 47rr 


Q in — 10 
• iO m 


suriace area 


V- — —r^ 


5 • 10-16 


cell volume 


A 


0.2 • 10-6 m 


thickness of extracellular space around cell 


Vextra — S l\. 


6 • 10-1"^ 


extracellular volume 


Vmax,l ~ 2 • 10 S 


6 • 10-19 kg s-i 


maximum speed of GLUT-mediated glu- 
cose transport 




1.2 • 10-19 kg s-i 


maximum rate of glucokinase activity 


'^max.22 


1.2 • 10-1^ kg s-i 


maximum rate of hexokinase activity 


Ki 


0.27024 kg m-3 


Michaelis-Menten constant of glucose 
transport 


K2 


1.8 kg m--^ 


Michaelis-Menten constant for glucoki- 
nase 


K22 


1.8 • 10-2 kg 


Michaelis-Menten constant for hexokinase 


Ka 


5.4 • 10-2 kg m-^ 


constant for the homeostatic loop control- 
ling glucokinase and hexokinase activity 


D 


7-10-10 m2 s-i 


diffusion constant of glucose in water 




0.9 kg m-^ 


glucose concentration in standard nourish- 
ing solutions 



Table 1 

Parameters used in the numerical tests with one isolated cell. Parameters are ex- 
plained and referenced in \^fb\ . except D [21] and A (see discussion in the main 
text). 

The nearly exponential transient is clearly visible in a plot of the mass flow 
into the extracellular space 



dm, 



extra 



dt 



-05 / ^extra 



inflow 



extra 



which is shown in flgure[T]for the initial conditions mj„ = 0.1 [Genv] Vin and 
f^extra = 0.1 [G env] Vextra- Becausc of this fast transient, many standard algo- 
rithms fail when integration is carried out on long time intervals with com- 
paratively long time steps: rather surprisingly, Mathematica itself fails badly 
(at least with the standard settings for NDSolve) even though it includes a 
stiffness detection procedure [19]. This failure becomes slowly apparent as 
the integration is carried out on longer and longer intervals; flgure |2]) shows 
the concentration [Gin] vs. time for a 1000 seconds interval and initial condi- 



8 



8.x 10- 



X, 6.x 10-13 

(50 

M 

I 4.x 10-13- 
1/1 



2.x 10 



-13 _ 







0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 

time (s) 

Fig. 1. Mass flow into the extracellular space (mass flow (kg s~^) vs. time (s)) 
obtained from the numerical integration of equations (|3]) and (|4| with Mathematica. 
A direct fit of this nearly exponential curve yields a decay time r ~ 62/is, very close 
to the estimate in the main text. 

tions rriin = 0.9 [Genv] Vi„ and rriextra = 0.1 [Genv] Vextra- instead of a smooth 
decay towards the equilibrium value the Mathematica solution shows some 
unexpected undulations. The companion plot for the concentration in the ex- 
tracellular space [Gextra\ (not showu) displays a large, unphysical peak close 
to the origin (this peak is about 100 times larger than the environmental con- 
centration). A similar integration performed on a longer interval (10000 s) 
returns an even worse result, as it produces negative and unphysical values 
for the concentration [Gj„] (see figure IS)). 



A careful analysis of the curve in figure [T] reveals that the fast transient is 
very close to an exponential, and therefore we expect the standard results 
on A-stabihty of integration methods [T6] to be appUcable to this and to 
similar nonlinear differential systems. In particular, we know that both the 
implicit Euler and the implicit trapezoidal method are unconditionally A- 
stable. Indeed the tests performed show that both methods converge even 
when the internal machinery of Mathematica - which uses a full array of 
sophisticated methods - is foiled by this stiff system. And although it is known 
that the accuracy of the implicit trapezoidal method (IT) is higher than that 
of the implicit Euler algorithm (IE) [IB] , in these tests IE fares better, since the 
IT method shoots below the equilibrium value of [Gj„] before swinging - very 
slowly back (see figure |4]), and displays unwanted oscillations of [Gextra\ that 
very gradually die out. The IT method also alternates steps above and below 
the equilibrium value of [Gextra], and these oscillations depend on the step size 
and resemble those observed in a closely related method, the Crank-Nicolson 
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■ 200 400 600 800 1000 

time (s) 

Fig. 2. Plot of the concentration [Gjn] (kg m^^) vs. time (s) obtained from the 
numerical integration of equations ^ and Q with Mathematica for a 1000 sec- 
onds long time interval. The gray horizontal line shows the equilibrium level. The 
internal stepping procedure (with the default values) fails and produces unexpected 
undulations. 
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Fig. 3. Plot of the concentration [Gm] (kg m~^) vs. time (s) obtained from the 
numerical integration of equations ^ and ^ with Mathematica for a 10000 seconds 
long time interval. The gray horizontal line shows the equilibrium level. The internal 
stepping procedure (with the default values) fails badly and leads to large, negative 
(unphysical) values of the concentration [Gin] ■ 
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Fig. 4. Plot of the concentration [Gin] (kg m"'^) vs. time (s) obtained from the nu- 
merical integration of equations ([S]) and Q with the implicit Euler method (dashed 
curve) and with the implicit trapezoidal method (dotted curve) for a 10000 seconds 
long time interval; in each case the step size is 100 s. The gray horizontal line shows 
the equilibrium level: the implicit Euler method reaches this equilibrium value with- 
out unwanted oscillations, while the implicit trapezoidal method shoots below the 
line and afterwards climbs very slowly back to the equilibrium value. We find that 
this climb depends on the actual stepsize, and it is faster for shorter time steps. 
Neither algorithm produces unphysical values (negative concentrations). 

method for partial differential equations [22|23] . 



5 Diffusion in cell clusters 



In the previous section we have considered just one cell, but as explained 
above, we wish to simulate large clusters of cells, and this means that we 
have to deal with large systems of coupled nonlinear differential equations. 
We begin with a very small cluster of just two cells, where glucose dynamics 
is described by the following equations for the masses 
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where the superscripts (A) and (i?) denote the two cells, and the parameters 
in the numerical tests are a straightforward extension of those used for a single 
cell and are listed in table 2. 



The geometric structure of the two cells is depicted in figure [5] this structure 
approximates the shape of two cells just after mitosis, and we distinguish six 
different regions 

• the two hemispherical cells; 

• two hemispherical extracellular spaces that are the interface between the 
cells and the environment; 
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psirameter 


value 






5 • 10-^ m 


cell radius 


A 


0.2 • 10-6 m 


thickness of extracellular space around cell 


Al = r 


5 • 10-6 m 


effective distance for calculation of flow 
into aiiQ out oi tne extraceiiuiar space ue- 
tween cells 


.9 = l-r-' 


3 ■ 10"^" m 


lolal surface area 


5(A) ^ 5(B) ^ 


1.5 • 10-1° m 


surface area of each hemisphere 


q(AB) _ „„2 


( .y • lu 111 


COIlXdCt (XTHdb OGtWGGIl Cello 




1^.1 n^io 

1.0 lU 111 


huliduc fxied ui cacii iiciiii&piici c 


y{A) _ ^(B) _ 2^ 3 
in in 3 


2.5 • 10-16 


cell volume 


y^Sra = vSL = 2-^A 


3 • 10-1^ m2 


volume of each hemispherical extracellular 

space 


^extra ^ 


1.6 • 10-1^ 


volume of extracellular space between cells 


max,l max,l 


3 • 10-19 kg s-i 


maximum speed of GLUT-mediated glu- 
cose transport between each cell and 
surrounding hemispherical extracellular 
space 


JAB) 
max,l 


1.6 • 10-19 kg s-i 


maximum speed of GLUT-mediated glu- 
cose transport between each cell and the 
extracellular space uetween ceiis 


^max,2 


1.2 • 10-19 kg s-i 


maximum rate of glucokinase activity 


Vmax,22 


1.2 • 10-18 kg s-i 


maximum rate of hexokinase activity 


Ki 


0.27024 kg m-3 


Michaelis-Menten constant of glucose 
transport 


K2 


1.8 kg m-^ 


Michaelis-Menten constant for glucoki- 
nase 


K22 


1.8 ■ 10"- kg m"'^ 


Micliaclis-Mciilcn couslaut. for licxt)kinasc 


Ka 


5.4 • 10-2 kg 


constant for the homeostatic loop control- 
ling glucokinase and hexokinase activity 


D 


7 • 10-1° s-i 


diffusion constant of glucose in water 


[Genv] 


0.9 kg m-^ 


glucose concentration in standard nourish- 
ing solutions 



Table 2 

Parameters used in the numerical tests with two cells. 
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Fig. 5. Schematic view of the geometric structure of the pair of daughter cehs 
described in the main text, where it is assumed that just after mitosis the cehs 
are (nearly) hemispherical and are surrounded by two hemispherical extracellular 
spaces - the interface with the environment - and are separated by one cylindrical 
extracellular space. Left panel: perspective view; right panel: cross section. For clar- 
ity the thickness of the extracellular spaces is exaggerated; the interface between 
extracellular spaces is actually negligible and is not counted in the equations (|9| to 

(lU). 

• the cylindrical extracellular space between the cells; this cylindrical space 
is in contact with the environment only through a narrow belt; 

• the (fixed) environment. 

This rough division of space contains the main elements of simulations where 
diffusion plays an important role, and later in this section we shall see how to 
expand it. 

We have carried out extensive numerical tests similar to those in the previous 
section, and essentially they confirm what we found with a single cell: 

• the Mathematica solution obtained with the default settings fails for long 
integration intervals; 

• the implicit Euler method and the implicit trapezoidal method both con- 
verge to the equilibrium value, but the implicit trapezoidal method displays 
strong and very undesirable oscillations. 

This small system with just two cells can be extended further, with the ad- 
dition of more cells and their extracellular spaces: this is equivalent to the 
discrete volume version of a diffusion problem (see, e.g., [23] for a review of 
volume discretization for the solution of diffusion problems), i.e., we must deal 
with a large system of equations (here we use concentrations p rather than 
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the masses because the ensuing analysis is shghtly easier to follow): 

dp, 
^ dt 



vj^ = Fa (p,; Pa) + DY, (p, - p,) g,, (14) 



where indexes denote both cells and extracellular spaces, Qab is term related 
to geometry (it is equivalent to the 5/ A ratio that we met previously), the 
Fq's describe facilitated transport (for both cells and extracellular spaces) and 
metabolic processes (cells only), and the diffusion term is actually included 
only in the case of extracellular spaces. The subscript (6) indicates that the 
sum is performed only on the set of cells h that are adjacent to cell a: since 
in a simulator like VBL the cells are in arbitrary positions in space [3j, and 
not on the usual cubic lattice, the number of neighbors is random (in a 3D 
configuration it fluctuates about an average of 12). 

We have already remarked in the previous section that the F^'s are slowly 
varying functions and that the stiffness of the differential system and the 
algorithmic stability problems are almost entirely due to the diffusion terms: 
for this reason we concentrate on a reduced version of the differential system 



(14) 



^. D ^^^^ 



dt Va 



It is well-known that differential systems like (15), obtained from the dis- 
cretization of diffusion problems on regular lattices, can be integrated by 
implicit methods like the Backward Differentiation Formulas (BDF) or the 
Crank-Nicolson algorithm (CN) [I7I22], and that in those cases both BDF 
and CN are unconditionally stable However in simulation programs like 
VBL, the equations are not discretized on a lattice, and may have a variable 
number of diffusion terms, as the sum runs over a random number of neighbors: 
are these implicit integration methods still stable in this more general setting? 
Fortunately the answer is yes, because it is possible to develop a variation of 
the standard Von Neumann stability analysis. The argument for BDF goes as 



follows: we start from the discrete-time version of ( 15 ) that corresponds to the 
BDF iterations 

•'a 



{b) 



where At is the time step, p^ = pa{nAt), and we take the test solution p" ~ 
A" exp {ik ■ Ta) that corresponds to the test solution used in the standard Von 
Neumann stabihty analysis, but without the assumption of a regular spatial 
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lattice. Substituting the test solution in the iteration formulas (16) we find 

, . DAt 



A = I + —r^ 2^ gab [l - cos (pab) 



Va 



{b) 



ab Sm (pab 



(17) 



where (j)ab = k- (rb — ra), and since J29ab{^ — cos (pab) > we see that |y4| < 1 

{b) 

and thus the BDF algorithm is unconditionally stable. 



The CN iteration formula is 
D 



=pl + At 



2Va 



{b) 



- Pa ) 9ab + 2^ [Pb 



Pa)9ab 



(b) 



and we can repeat the steps followed for the BDF algorithm and find 



A 



1 - 



DAt 

Va 



E{6> 9ab{^- COS (j)ab) +i^Y.{b) 9ab siu (pab 



1 + 



E(b> 9abO- - COS. 



>ab 



;DAt 

' Va 



E(6> 9ab sin (pab 



(19) 



and just as before we see that < 1 and that the CN algorithm is uncondi- 
tionally stable. 

To test the BDF and the CN algorithm in a realistic biophysical setting, we 
have extended the two-cell model to a ID string of cells, where we assume 
that each cell is surrounded by its own extracellular space, and that adjacent 
extracellular spaces can exchange material by diffusion. This means that the 
diffusion problem becomes slightly more complex (although the previous sta- 
bility analysis still holds) and the equations for each cell/extracellular space 
are 



Vc^ = M{pc)+T{p,,Pc) 

vj^ = Dj2{pb- pc) 9bc - T {pc, pc) 
"'^ (6) 



(20) 
(21) 



where 



2 

Mg {PG,C) ■ 



{K2 + pG,c) {Ka + pG,c) 
2 

VimiX,22PG,C ^22) 



(^22 + PG,C) (Ka + pG; 



C) 



Tg {pG,C, pG,c) = j^— j^— 23 

Ki + Pg,c Ki + pg,g 
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Fig. 6. Contour /density plots that show different aspects of the solutions of equa- 
tions (20 21 ) for a string of 100 cells, that spans a spatial range (— 500//m, 500/um). 



Space runs along the x-axis and time runs upwards along the y-axis with timesteps 
At = 2 s, and concentration values are also mapped on a gray scale (black corre- 
spond to the minimum value, white to maximum). The panels show: a) the contour 
plot of the glucose concentration in the extracellular spaces obtained with the BDF 
method; b) contour plot of the glucose concentrations inside cells, for cells in the 
reduced spatial range (— 500/im, — 200/im), for t = s until t = 60 s, obtained with 
the BDF method; c) same as b, but in this case the numerical solution has been 
calculated with the CN method: notice that the contour lines are appreciably dis- 
torted close to the boundary; d) contour plot of the glucose concentrations in the 
extracellular spaces in the reduced spatial range (— 500^m, — 200/im), for t = s 
until t = 60 s, obtained with the CN method: notice the banded structure that is 
due to the oscillatory instability of the CN method. 

are the functions that describe the internal metabolic activity (conversion of 
glucose to G6P) and the transport processes into and out of cells and that 
we met earlier both in the one-cell and in the two-cell models. The upper- 
case subscripts denote cells and the lowercase subscripts denote extracellular 
spaces, and all the model parameters have already been listed in table 2. The 
space coordinate (cell positions) span the range {xmin,Xmax), with the bound- 
ary conditions u{xmin,t) = u{xmax,t) = constant: thus this is a model with 
both diffusion of glucose and local absorption and conversion into G6P. 

The BDF and the CN algorithm are equivalent to the implicit Euler and to the 
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200 400 600 800 1000 



time (s) 

Fig. 7. Plot of the glucose concentration in the leftmost extracellular space in the 
string of cells described in the text, obtained numerically with the CN method for 
the same problem of figure [6j In this solution there is a marked high frequency 
spurious oscillation that depends on the time step and that very slowly fades away. 




200 400 600 800 1000 



time (s) 

Fig. 8. Plot of G6P concentration at two positions (close to the boundary and in 
the middle of the string of cells), obtained numerically with the BDF method for 
the same problem of figure [6| using the differential equation for G6P in [3] , solved 
with the implicit Euler method; the concentration of G6P is almost independent of 
position along the string and the two curves are indistinguishable on this scale. 
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Fig. 9. Plot of G6P concentration along the string of cells at t = 1000 s, obtained 
numerically with the BDF method for the same problem of figure [6j using the 
differential equation for G6P in jl], solved with the implicit Euler method. The 
concentration is lowest in the center, and it mirrors a similar curve for glucose 
concentration. 

implicit trapezoidal method, respectively (they are actually the extensions of 
these algorithms to partial differential equations), and thus - unsurprisingly 
- they share with their counterparts both advantages and disadvantages. In 
particular in our tests on this particular diffusion problem we find that: 

• both BDF and CN eventually converge to the same equilibrium solution; 

• the BDF method requires a larger number of function evaluations to reach 
the required accuracy; 

• the CN method displays unwanted oscillations, similar to those found in the 
case of the implicit trapezoidal method. 

Figure |6] shows contour plots obtained with the BDF and the CN methods, 
and figure [7] shows the CN solution close to the boundary vs. time: while the 
global behavior of the solutions is the same, the slow-dying oscillations of the 
CN method close to the boundary are clearly visible. 

The oscillatory behavior of the CN method is due to the reduced damping 



of the high-frequency spatial modes (see the formula for A, eq. (19)), and 
it is known to sometime spoil the quality of the CN solutions |23j. These 
numerical tests indicate that the BDF algorithm - as an extension of the 
implicit Euler algorithm - is a good choice for an integrator in large scale 
biophysical simulations. 
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Fig. 10. Simulated tumor spheroid in VBL. The simulation started from a single 
cell and this cluster of cells corresponds to more than 6 days of simulated time. 

6 Conclusions 



We are developing a simulation program, VBL (Virtual Biophysics Lab) where 
we aim to include a basic description of biochemical and biomechanical fea- 
tures, to simulate cell clusters, and that should eventually be a numerical 
model of tumor spheroids, a useful and important in vitro model of solid tu- 
mors [25]: the stakes are high and this is just one of several attempts to use 
mathematical and physical methods to understand the biochemical and me- 
chanical properties of tumor spheroids [26]l?ril28]l2ro]l8T]^^ 



In our modeling effort we proceed in a partly phenomenological way that leads 
to simple parameterizations: in exchange, we achieve a huge reduction in com- 
putational complexity and a considerable reduction of the space-time scale 
problems that affect simulations aimed at calculating the properties of macro- 
scopic objects starting from microscopic models. We are in an advanced phase 
of development of the program, and we have already included cell metabolism, 
growth and proliferation and the extracellular environment. In the previous 
sections we have discussed the stability properties of the integrators in VBL 
and in similar programs: we have performed tests on a simple model system 
and some of the results are interesting on their own, like, e.g., the time evo- 
lution and the distribution of G6P in a ID string of cells (see figures 10 and 



20 



60 



40 



20 







-20 



-40 



-60L 



6 22: 56: 40 



-60 




-40 



-20 



20 



40 



60 



Fig. 11. Distribution of dead cells (dark gray) inside the simulated tumor spheroid 



of figure 10 Here cells are represented as semitransparent balls, and dead cells in 
the core become visible. The numbers on the axes are /um. 



11). The 3D part of the program, i.e. geometry and biomechanical interac- 
tions, is also included in the prototype version of the program, as well as a 
description of the extracellular microenvironment. We can already simulate 
large populations of dispersed cells, like those in the culture wells used for in 
vitro growth, and we have produced numerical estimates that are in excellent 
qualitative agreement, and in good quantitative agreement, with experimen- 
tal data |^|f5|l6] . Figure 10 shows one simulated spheroid and figure 11 is an 



alternate display that shows the distribution of dead cells inside the same 
spheroid. 



In the present stage of development of the program we are cleaning up the 
biochemical part of the simulation and building a close integration with the 
diffusion algorithm to achieve a robust, stable program. Because of discrete - 
and sometimes random - events in the cells' lives, the evolution is only partly 
described by differential equations. Moreover - because of cell proliferation - 
there is a variable number of equations. On the basis of the considerations 
and the numerical tests described in this paper we have decided to settle on 
the implicit Euler algorithm and its extension, the Backward Differentiation 



Formula. The simulations that produced figures 10 and 11 were based on much 
simpler (and unstable) explicit Euler integration steps, and stability problems 
showed up at an early stage: soon we shall be able to carry out higher quality 
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simulations based on the much more stable and reliable methods described 
here. 

Eventually we may have to tackle efficiency and speed as the simulations be- 
come larger and more time-consuming. However the choice of BDF is open to 
improvements in speed, since - as we have seen in the last section - biochem- 
ical changes are comparatively slow, and stiffness is brought about mostly by 
the diffusion terms in the equations: this means that we may eventually use 
implicit-explicit methods like those described in |13] and this may help us re- 
tain the advantages of the implicit BDF algorithm and greatly increase speed 
- but this is something that we shall study in a future work. 
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